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Level set methods are versatile and extensible techniques for general front 
tracking problems, including the practically important problem of predicting the 
advance of a firefront across expanses of surface vegetation. Given a rule, em- 
pirical or otherwise, to specify the rate of advance of an infinitesimal segment 
of firefront arc normal to itself (i.e., given the firespread rate as a function 
of known local parameters relating to topography, vegetation, and meteorology), 
level set methods harness the well developed mathematical machinery of hyper- 
bolic conservation laws on Eulerian grids to evolve the position of the front in 
time. Topological challenges associated with the swallowing of islands and the 
merger of fronts are tractable. 

The principal goals of this paper are to: collect key results from the two 
largely distinct scientific literatures of level sets and firespread; demonstrate 
the practical value of level set methods to wildland fire modeling through nu- 
merical experiments; probe and address current limitations; and propose future 
directions in the simulation of, and the development of decision- aiding tools to 
assess countermeasure options for, wildland fires. In addition, we introduce 
a freely available two-dimensional level set code used to produce the numerical 
results of this paper and designed to be extensible to more complicated configu- 
rations. 

Keywords: Wildland firespread, level set methods, Multivac software 



1 Introduction 



Wildland fire modeling has received attention for decades, due to the sometimes 
disastrous consequences of large fires, and the trer nendous cost s of often ineffec- 
tual, possibly even counterproductive firefighting Pvnd . 2004 1. For the practi- 
cally important scenario of wind-aided firespread, one seeks a computationally 
efficient model, useful not only offline (for pre-crisis planning, e.g., placement of 
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access roads, firebreaks, and reservoirs, and scoping of fuel-reduction burning, 
and post-crisis review, e.g., personnel training, litigation), but also during a cri- 
sis (i.e., real-time guidance for evacuation and firefighting) . For computational 
efficiency, such that the benefits of ensemble forecasting [Palmer et al. . 2005l | 
are readily accessible from a model, advantage should be taken of the inherent 
scale separation of: (1) the kilometer-and-larger, landscape-dominated scales 
of the local atmospheric dynamics; and (2) the one-meter-and-smaller scales of 
the local combustion dynamics. Even with advanced techniques and access to 
exceptional contemporary computing facilities, numerical simulations (of turbu- 
lent flows) that proceed from fundamental principles are challenged to resolve 
accurately in real time phenomena with spatial scales spanning much more than 
two orders of magnitude [H.R. Baum, private communication]. Thus, the feasi- 
bility of a direct numerica l simulation encompassing the m ultivaried processes o f 
wildland fire propagation [Coed , may be decades off [Jenkins et al.l . l200lj . 

Moreover, at least many attempts (albeit usually problematic) at parameter- 
ization of subgridscale phenomena in terms of gridscale variables have been 
undertaken by meteorologists for cumulus convection, turbulent transport, and 
radiative transfer. However, meteorologists have extremely limited experience 
with the parameterization of combustion dynamics for weather-dependent wild- 
land firespread; even if such parameterization be possible, it remains unknown. 
Furthermore, data collection in wildland fires is so piecemeal, irregular, and of 
uncertain accuracy that, for many years to come, the data better suit reinitial- 
ization of a simplistic model than assimilation into an ongoing calculation with 
a highly detailed model. 

Accordingly, in this study, attention is focused on a minimalist treatment 
of the firefront, idealized as an interface between expanses of burned and un- 
burned vegetation. This treatment is consistent with the typically limited, only 
gross characterization available for the vegetation at issue, since the vagaries 
of ignition events are difficult to anticipate, and maintaining an updated in- 
ventory for the huge area of wildlands in (say) the USA is daunting. This 
simplistic interfacial approach to the fire dynamics, easily executed in minutes 
on a laptop given the requisite meteorological and other input fields, reserves 
computational resources for the difficult, more critical, and mostly yet-to-be- 
undertaken landscape-scale weather forecasting targeted for real-time wildfire 
applications. 

The upshot is that simple persistence models are adopted for the wind field 
(and thermodynamic variables) in the study undertaken here. Also, attention 
is limited to a one-way interaction between the meteorology and the firespread, 
though future extension to two-way interaction by use of an iterative procedure 
may be envisioned. Simplistic modeling still may provide the key macroscopic 
fire behavior sufficiently accurately for practical purposes (including estimates 
of smoke and pollutant generation), even for circumstances for which the sim- 
plification is not formally justifiable. In fact, observational data of wind-aided 
firefront progression in wildland are today typically sparse, so that not much 
more than the output of a simplistic model can be meaningfully validated and 
tuned. Moreover, the use of relevant mathematical methods to perform model 
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selection, to carry out efBcient parameter estimation, and to account for the 
uncertainty in predictions is facilitated by focusing on less detailed models with 
fewer parameters. In this paper, we mainly address the first step, which is to 
achieve proper forward simulations. 

One of the most widely used models was devised by Rothermel Rothermei 

to predict the rate of firespread, with focus on the head of a wind-aided 
fire. Because predictions of the Rothermel treatment have been found to be at 
odds with some observations, efforts to improve this spatially one-dimensional 
semi-empirical treatment, and to supplement t he data upon w hich it is based, 
have been undertaken, especially in recent years Carltonl . l2003l |. Extension from 
a focus exclusively on the head of the fire seeks to evolve the configuration of 
the entire fire perimeter, possibly of multiple fire perimeters. In this study, and 
typically, the firefront, even a moderate fraction of an hour after a localized 
ignition in fire-prone vegetation, is taken to be a closed curve projected on a 
plane (the ground may n ot be flat). Such simulations of firespread have been 
performed Finney . 1998| with the so-called marker technique, which discretizes 
a front into a set of marker particles, and advances the front through updates of 
the particle positions. Parenthetically, as a problematic step, the updating by 
Finney takes each marker on the front to evolve identically to an idealization 
of how a front evolves from a single isolated ignition site in an unbounded 
expanse of vegetation, in the presence of a wind. In any case, ev en though 
applied projects have supported software development Finney, 1998| . still from 
a computational point of view, only a f ew, largely equival ent methodological 

In t his paper, we 
1999t to calculate 



ecu 

developments have been undertaken [e.g. lAndre et al.L 1200' 



Sethia: 



apply level set methods Osher and Sethianl . Il988l : 
firefront evolution. 

In Section [21 we introduce wildland fires pread models, especially a semi- 
empirical, equilibrium-type model proposed in F endell and Wolfn 200 1| for wind- 
aided firespread across surface-layer, chaparral-type, burning-prone vegetation. 
(In commonly adopted equilibrium-type models, the firespread rate depends on 
only the parametric values holding locally and instantaneously, so the firespread 
rate is taken to adjust indefinitely rapidly to any temporal and spatial change.) 
Section [3] provides a brief introduction to level set methods. Section 2] describes 
the Multivac level set package that has been applied in this paper to the fire- 
spread problem. A quick description of its performance is presented in Section[5l 
Finally, results of firespread simulations with different idealized environmental 
conditions are reported in Section [6l 



2 Front Propagation Functions for Wildland Fires 

Even if theory and/or measurement furnished complete, perfect knowledge of 
the topography, vegetation, and meteorology at a site at a given time (e.g., fur- 
nished the locally pertinent values of all parameters in functional forms capable 
of representing these three types of input), still one currently possesses very 
incomplete, imperfect knowledge of the "rules" that would yield the physically 
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observed rate of firespread from the input. Achieving knowledge of firespread 
"rules" sufficiently accurate for practical purposes may well lag emplacing means 
for observing and collecting exhaustive input data. 



As already noted, a fire-growth simulation such as FARSITE [Finnevl . Il998l | 
seems unlikely to reach its potential as long as it seeks to describe the rate of 
firespread at all orientations to the direction of the sustained low-level ambi- 
ent wind from sp read-rate modeling focused on the direction of the wind [e.g. 



Rothermell |1972| | . On the other hand, posing a different rule for the spread rate 



at every possible orientation to the wind defeats the goal of simplicity. 
2.1 Wind-aided wildland fire spread 

Fendell and Wolff Fendell and Wolfi . 2001*1 addressed this dilemma in develop- 



ing a model dedicated to wind-aided wildland fires that spread rapidly over level 
terrain with dry, moderately sparse fuel, taken here to be uniformly distributed 
to permit concentration on wind effects. Parenthetically, for consistency with 
modeling in which the firefront is idealized as an interface moving according to 
a semi-empirical rule, only a minimal amount of information about the surface- 
layer fuel is required, mainly the mass loading consumed with firefront passage 
( "available" -fuel loading). 

The Fendell and Wolff model focuses on front velocities at the rear of the 
front (where propagation is against the wind), at the head of the front (where 
propagation is with the wind), and on the flanks (where propagation is across 
the wind direction) - see Figure [T] The firespread velocities primarily depend 
on the wind velocity U . At the rear, the front advances relatively slowly against 
the oncoming wind, since hot combustion products tend to be blown over an 
already burned area. The velocity at the rear is denoted e{U). At the head, 
the velocity h{U) is relatively large, since hot combustion products tend to 
be blown over a yet-to-burn area, in which discrete fuel elements are heated 
toward ignition by convective-conductive transfer. Both analytic modeling and 
laboratory experim ents have shown that h{U) is roughly proportional to ^/U 



Wolff et al.l . Il99l| . At the flanks, the (spread-aiding) wind component along 
the normal to the front is zero, but observationally the front advances faster 
than in the absence of wind. As a speculation, a more meticulous treatment 
would find that, at the nominal flank, the configuration is convoluted, and 
firespread is alternately with and against the wind. Of course, were the wind 
direction constant, limiting attention to the head would seem adequate, but, 
in fact, change in wind direction may (rapidly) result in an interchange of the 
locations of the flank and head - an interchange sometimes associated with 
tragic consequence for firefighters. 

The velocities (the termi nology henceforth adopte d, for brevity, in place of 



firespread rates) proposed in iFendell and Wolfj [200l| are 



e(C/) =eoexp(-eiC/), /([/)= + ciC/exp(-c2C/), hiU) = eo + aVU, (1) 
where Eq, Ei, ci, C2 and a are parameters (with readily inferred dimensionality) 
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Figure 1: Fendell and Wolff model introduces velocities at t he rear (against the 
wind), at the head (in the wind direction) and at the flanks. 
200 1| 



Fendell and Wolfi. 



depending on the mass loading of fuel and other parameters characterizing the 
fuel bed, but independent of U. 

The velocity is then provided at any point on the front through a "trigono- 
metric interpolation" : 



F{U,9) = fiU siire)+h{U cos" 9) if \0\ < f, 
F{U, 0)=f {U sin™ 9) + e{U cos^ 9) if 16*1 > f , 



(2) 



where 9 is the angle between the wind direction and the normal to the front. 
We set m — 2. In this paper, parameter n is set to | and is significant since it 
determines the overall shape of the front from the flanks to the head. 
To summarize, the velocity is, for all U £ K+ and 6* s] — tt, 7r[, 

F{U,0) = eosm^ 9 + cqU am^ 9 ex:p {-^ciU am^ 9) 

Eq cos"^ 9 + aVu cos"- 9 if |0| < f (3) 

So cos^ 9 exp (-eiJ7 cos^ 9) if |6l| > f ' 

2.2 Simplified model 

Based on the numerical experiments carried out with the level set co de Multi- 



vac (Section |4|), the model ([3]) proposed in IFendell and Wolfj [200l| has been 
modified. First, the parameter n has been set to | instead of 1. Second, the 
model has been simplified without losing its main features, primarily the overall 
shape of the firefront. The new model reads 

FiU,9)=ea + ciVUcos"9 if 1^1 < f, 



F{U,9)^ea{a + {l-a)\sm9\) H \9\ > 



2 ' 
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where a £ [0, 1] is the ratio between the velocity at the rear (a£o) and the 
velocity at the flanks (eo). Velocities at the rear and at the flanks no longer 
depend on the wind, since their dependence on the wind speed is hard to model 
accurately and has little impact on the overall front location. The velocity at 
the head is the same as in the "full" model ([3]). 

The simplified model is easier to tune, either via direct trials or with sys- 
tematic methods for parameter estimation (which may require derivatives of 
the model with respect to its parameters). All results in this paper are for the 
simplified model. However, results for the "full" model would appear roughly 
the same. 

3 Level Set and Fast Marching Methods 

First introduced in lOsher and SethianI |l988l | , level set methods are Eulerian 
schemes for tracking fronts propagating according to a given speed function. In 
this section, we explain basic features of the level set methods used for firespread 
modeling. 

3.1 Mathematical basis and technique 

3.1.1 Definitions 

Assume the front evolves from the initial time i = to the final time t = Tf. 
For all t e [0, Tf], the front at time t is the set of points (in E^) T{t). We define 
To = r(0) as the initial front. 

For all t G [0, Tf], each point X G T{t) with a well-defined normal moves in 
the direction normal to the front with a given speed F{X,r,t). Notice that F 
may depend on the position, on the time and on local properties of the front 
itself (certainly the normal direction, not always defined, and possibly the local 
curvature or other properties). 

The problem is to approximate T : [0, Tf] K^, given Fo and F. 

3.1.2 Strategy 

The main idea is to evolve a function ip : x [0, Tf] R such that 

yte[o,Tf] r(t) = i^x eR^ /ip{x,t) (5) 

(fi is called the level set function. At any time, the zero level set of (p is the 
front itself. A priori, ip could be any function satisfying equation ([5]). However, 
some assumptions (e.g., smoothness) and practical issues (e.g., initialization of 
(f) make it convenient to define (p as the signed distance to the front. 

Then, if d is the Euclidean distance on R^, we define, for any given curve 
T, the distance to T: 

VxeR^ dr{x) =mmi^d{x,P)/p (6) 
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Hence the signed distance ip for all x £ and t £ [0, T/]: 

I ) d.r(t){x) if X lies outside the front r(t) , , 

" \ -drit){x) if X lies inside the front r(i) ' ^ ' 

It can be shown that ip obeys the equation 

VxeR^ \ft£[0,Tf] ^t{x,t)+F{x,^{-,t),t)\\\/Mx,t)\\^^0, (8) 

where the velocity F is now defined everywhere in an d depend s on th e front 
through its dependence upon p. Details may be found in Sethian 1999| . 



Recall that p{-, 0) is known as well as To; ip{0) is the signed distance to F, 
pjv / — J '^^oix) if X lies outside the front Tq 



0- 



^ ' ' [ ~dro{x) if a; lies inside the front Fq ' 

Equations ^ and ([5]) define the initial-value problem that is to be solved. 
Zero level sets of ip yield the front points. 

This nonstationary problem involves the Hamilton- Jacobi equation ([5]) . There 
may be multiple solutions to this equation. P.-L. Lions and M. G. C randall 
defined the so-called "vis c osity solution" of Hamilton- Jacobi equations [LionsL 
1982t ICrandall and Lioil^ . Il983l |. which turns out to be the unique physical so- 



lution for which we search. Under given assumptions (mainly on the speed func- 
tion F), existence and uniqueness of the viscosity solution of the problem ([S])-® 
can be proved. 



3.2 Advantages and disadvantages of level set methods 

Several methods may be relevant to simulate the propagation of fircfronts. One 
may want to use marker techniques, in which the front is discrctized by a set 
of points. At each time step, each point is advanced according to the speed 
function. This Lagrangian methodology leads to low-cost computations, but 
requires care in the handling of topological changes. 

Volume-of-fluid methods represent the front by the amount of each grid-cell 
that is inside the front. In each cell, the front is approximated by a straight 
line (horizontal or vertical, in most methods). Such methods can deal with 
topological changes, but the front representation can be inaccurate. In wildland 
firespread, the direction normal to the front is crucial because of the wind- 
direction-dependent speed function (see Section [2]). 

Level set methods automatically deal with topological changes that occur in 
wildland firespread, such as fronts merging and front convergence (in connec- 
tion with unburnt "islands"). The level set description enables a fair estimate 
of the normal to the front, making it well suited to the fire propagation problem. 

However, level set methods have disavantages. First, they embed the front 
in a higher-dimensional space. Helpfully, the narrow band level set method 
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Adalsteinsson and Sethian . 1995| is an efficient algorithm which almost de- 
creases the problem d imension by o ne. Moreover, when it can be used, the 
fast marching method Sethian . 1996l | provides a highly efficient algorithm. 

The main reservation may be the lack of proof of convergence of numerical 
schemes for certain problems. For a given class of speed functions, th e prob - 
lem dH])-® may routinely be solved numerically Crandall and Lion^ . 1984 1. 
However, no proof of convergence in mesh parameter or time step is yet avail- 
able for some situations. 



3.3 Quick review of numerical approximations 

Numerical approximation to solutions of Hamilton-Jacobi equations is closely 
related to numerical approximation to hyperbolic conservation law^. The point 
is to introduce a numerical Hamiltonian to approximate the Hamiltonian H = 

Crandall and Lions have proven that, for given Hamiltonians and initial con- 
ditions, a consistent, monotonic and locally Lipschitzian numerical Hamiltonian 
yields a sol ution that converges to th e vis cosity solut i on. F ormal results may 
be found in Crandall and Lion^ [l984 1 and ISouganidi^ 1985j . 



In one dimension, ipt+H(\/xif) = may lead to the following approximation: 
.-«..--A<.(^.^). (10, 

For instance, if the Hamiltonian is not convex, the Lax-Friedrichs scheme 
may be used; then, the numerical Hamiltonian is 

Va,6eR g{a,b)=H(^)-^^—^, (11) 



where the monotonicity is satisfied on [—R,R] if = max \H'{a)\. 

—R<a<R 

Several schemes have been developed, from simple and efficient schemes as 



that o f Engquist-Osher to high-order essentially nonoscillatory schemes |Osher and Shu 



1991 



3.4 Overview of complexity issues 

Let the mesh (in R^) be orthogonal with M points along each direction. Assume 
that the front is described by 0{AI^~^) points. The narrow band level set 
method makes it sufficient to update the level set function only in a narrow 
band (of width k) around the front. For each time step, the complexity of the 
algorithm is therefore 0(A:Af^~^). 

For an explicit temporal discretization the number of iterations is related to 
the Courant-Friedrichs-Lewy condition. Along x, the Courant number must be 

^Notice that, from equation (fix satisfies a hyperbolic conservation law in the one- 
dimensional case. 
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less than 1: 

< 1- (12) 

Usually, controlling the accuracy of approximation is subordinate to space 
discretization, which means that the time step is adjusted so that the Courant 
number is taken close to 1. 



Calculations may sometimes be sped up by reformulating the level set prob- 
lem as a statio nary problem. This leads to the so-called fast marching method 
Sethian . 199^ . Nevertheless, restrictions on the Hamiltonian prevent the use 
of this technique for some appl ications. The work of Sethian a nd Vladimirsky 
has overcome some limitations Sethian and Vladimirskvl . 12001 1. but restrictive 
conditions still remain (e.g., convexity of the Hamiltonian). 



4 Code 

4.1 Introduction to the Multivac level set package 

Multivac is a level set package freely available (under the GNU GPL license) 
at 'http : / / vivienmallet .net/fronts/. It is designed to be both efhcient and 
extensible, so that it may be used for a large range of applications. To achieve 
these goals, Multivac is built as a fully object-oriented library in C++. 

Multivac was designed independently of the firespread application described 
herein, but easily enabled firespread simulations, and is presently distributed 
with firespread-motivated fu nctions. It has also be en used in modeling the 
growth of Si-based nanofilms Phan and Mallet , 2003l | and image segmentation. 



The latest stable version available at the time of submission is Multivac 1.10. 



4.2 Structure 

The modularity of Multivac comes from its object-oriented framework, in which 
the main components of a simulation have been split into an equal number of 
objects. A simulation is defined by the following objects: 

• the mesh] 

• the level set function; 

• the velocity, which provides the propagation rate of the front according to 
its position, its normal, its curvature, and the time; 

• the initial front; 

• the initializer, which manages first initializations and initializations re- 
quired by level set methods (e.g., the narrow band reconstruction); 

• the numerical scheme, which advances the front in time; 
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• the output management. 

For each item, a set of classed with a common interface is available. For 
instance, several speed (i.e., propagation rate) functions are available through 
several classes, e.g. CConstantSpeed or CFireModel. All speed functions have 
the same interface, which allows users to define their own speed function on 
the same basis. The user principally provides speed rates as a function of the 
position, the time, the normal to the front and the curvature (these values are 
computed by Multivac itself). 



4.3 Calling sequence 

The whole is managed by an object of the class CSimulator. This object simply 
calls the initializer to perform the first initializations. Then it manages the loop 
in time (or iterations, in the case of the fast marching method) into which the 
numerical scheme is called to advance the front. The initializer is called again 
to reinitialize the signed distance function for the new step, and the object 
dedicated to post-processing requirements is called to save any needed data. 

In each step, objects communicate with one another through methods (i.e., 
functions) of their interface. For example, the velocity object provides speed 
rates to the numerical scheme. 



4.4 Overview of available classes 

Multivac package (version 1.10) includes several classes which are listed in Ta- 
ble [H 



4.5 Other strengths, limitations and future work 

Multivac takes advantage of C++ exceptions to track errors, and several de- 
bugging levels are defined, from a safe mode, in which all is checked, to a fast 
mode, in which performance is the primary concern. 

There are currently two main limitations. First, Multivac deals only with 
uniform orthogonal meshes. H owever, extensions of leve l set methods to un- 
structured meshes exist (e.g., Barth and Sethiaiil [l998|) and they could be 



implemented within the Multivac framework. Adaptively refined meshes are 
also accommodated with additional mathematical complexity, though the im- 
plementation effort would be substantial. Second, Multivac deals only with 
two-dimensional problems. 

Work is planned to allow inverse modeling (parameter estimation based on 
data assimilation) within the framework of Multivac. The main idea is to replace 
the class CSimulator with a class dedicated to inverse modeling. Preliminary 
results show the framework extendibility, but this capability is not yet available 
in distributed versions. Future versions should include this feature, based on an 
innovative method for integrating sensitivities along with the front itself. 



class is a user-defined type, in the manner of structures in C. Classes encapsulate data 
(called attributes) and functions (called methods). 
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Category 


Available classes 


Mesh 


Orthogonal mesh 


Level set function 


Defined on an orthogonal mesh 


Velocity 


Constant speed 
Piecewise constant speed 
Fire model 
Simplified fire model 
Image intensity 
Image gradient 


Initial front 


Circle 
Two or three circles 
One or two circles with an island inside 
Front defined by any set of points 


Initializer 


Basic initialization (no velocity extension) 
Extends the velocity with the closest 
neighbor on the front 


Numerical scheme 
(narrow band) 


Engquist-Osher, first order 
Lax-Friedrichs, first order 
Engquist-Osher, ENO, second order 
Chan-Vese algorithm [Chan and Vese. 2001] 


Numerical scheme 
(fast marching) 


Engquist-Osher, first order 



Table 1: Basic classes available in Multivac 1.10. 
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5 Complexity and Convergence Studies 



5.1 Convergence studies 



In this section, we report convergence studies th at are necessary to vahdate 
the code. As in lAdalsteinsson and SethianI [1998|, tests are carried out for a 
circle that expands in time with a unitary velocity. Details of the simulation 
are summarized in Table [H 



Data 


Value 


Comment 


Domain 


n = [0,3] X [0,3] 




Initial front 


Circle 




Circle center 


= (1.5,1.5) 


Domain center 


Initial circle radius 


^initial — 0.^ 




Final circle radius 


r final = 0.9 




Velocity 


F = 1.0 


Constant 


Duration 


Tf = 0.4 




Time step 


At = 10-4 





Table 2: Simulation test-case. 
We introduce three norms. The first is 

^spatial ~ simulated ~ ffinal\) 

where r simulated is the simulated radius, estimated as follows: 

1 



^simulated 



card(rrf) 



X! d((a;,y), (xc,yc)) 



(13) 
(14) 



where Yd is the discretized front as returned by the simulation (at time Tf) and 
d is the Euclidian distance. 

Additionally, if Ttrue(x, y) is the time at which the front is supposed to reach 
the point (x, y): 



p2 



\ 



1 



card(rd) 



{Tf -Ttrue{x,y)f 



The last norm is an infinity norm: 

^Ume ^ , max |T/ - Ttrue{x, 2/) I 



(15) 



(16) 



Table [3] shows results for the first-order Engquist-Osher scheme with the 
narrow band method. The width of the band is 12 cells and the front lies within 
a band whose width is 6 cells. 

The first-order Lax-Friedrichs scheme and the second ENO Engquist-Oshcr 
scheme were also checked successfully. As for the second-order scheme, the full- 
matrix method, that is, without the narrow-band restriction, was used because 
the front reconstruction destroys the second-order accuracy. 
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Ax = Ay 




^Ipatial (xlO^) 


eiime (xlO^) 


etT^e (xlO^) 


0.01 


301 


1.634 


1.753 


2.377 


0.005 


601 


0.855 


0.901 


1.191 


0.0025 


1,201 


0.460 


0.474 


0.600 


0.00125 


2,401 


0.244 


0.247 


0.299 



Table 3: Errors versus spatial discretization. 



5.2 Complexity issues 

Multivac was compiled under Linux with GNU/gH — h 3.3, and the reference 
simulation (see Table [2|) was launched on a Pentium 4 running at 2.6 Ghz. The 
width of the narrow band was 12 cells and the width of the inner band, in which 
the front lies, was 6 cells. If = Ny — 1001 (one million cells), the 4,000 
iterations were achieved in 14 s. 

The complexity of the narrow band level set method is close to 0{N), where 
N = Nx = Ny . Table d] shows that linear complexity of the method is not 
observed. Instead, the complexity seems to be 0{N'^). This is the complexity 
of the suboptimal algorithm currently used to rebuild the front. Moreover, the 
number of front reconstructions increases with the mesh refinement since the 
width of the narrow band does not change. 



Aa; = Ay 




Timings (s) 


0.03 


101 


0.4 


0.015 


201 


0.9 


0.01 


301 


1.6 


0.0075 


401 


2.6 


0.006 


501 


4.0 


0.005 


601 


5.6 


0.004285714 


701 


7.4 


0.00375 


801 


9.5 


0.003333333 


901 


11.9 


0.003 


1001 


14.1 



Table 4: Timings versus spatial discretization. 

6 Applying Level Set Methods to Firespread Ap- 
plications 

6.1 Method and numerical scheme 

The speed function ^ introduced in the level set equation ^ provides an 
Hamiltonian with nontrivial dependencies. Because of these dependencies (par- 
ticularly the non-convexity of the Hamiltonian) , neither the fast marching method 
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nor its extension to anisotropic problems can be applied. The narrow-band level 
set method is more relevant. 

A highly accurate numerical scheme is not required for the investigations 
reported here. The discrepancies between the numerical simulation and the 
exact solution should be considered in the context of other approximations: the 
model itself is simplistic; input parameters such as wind speed or fuel density 
are typically not accurately estimated; the location of the initial front introduces 
further uncertainties. A first-order scheme suffices for our purposes. 

Since the Hamiltonian involved is not convex with respect to spatial deriva- 
tives of the level set function, the first-order Lax-Friedrichs scheme (refer to 
equation (jlip ) is well suited. To minimize introduction of diffusivity, a local 
Lax-Fricdrichs scheme may be used as well. 

As previously advocated, the timestep Ai is chosen according to the Courant- 
Friedrichs-Lewy condition (fT2|) : 

01/S.x 

where a < 1; a is not kept constant in the tests that we undertake. Nevertheless, 
the Courant-Friedrichs-Lewy condition is estimated at every iteration with an 
(a priori) approximation to max|_ff'| along x and y, which leads to: 

At < (18) 



a(m + l)\/U 

The main characteristics of the simulation, including model parameters (refer 
to equation ([3])), are gathered in Table [S] 



Parameter 


Value 


n 


1.5 


U 


100 


a 


0.5 


eo 


0.2 


a 


0.5 


Ax 


3 • 10"^ 


Ay 


3 • 10-3 


At 


10--* 


Tf 


0.1 



Parameter 


Value 


Domain 


9. = [0,3] X [0,3] 


Initial front 


Circle 


Circle center 


(1.5,1.0) 


Initial circle radius 


^initial = 0.5 


Velocity 


F = 1.0 


Duration 


Tf = 0.1 


Time step 


A< = 5 • 10-5 


Spatial discretization 


N^ = Ny^ 1001 



Table 5: Parameters and their default values. 



6.2 Results 

The simulation described by Table [5] is shown in Figure [2j The figure shows 
snapshots of the front, initially circular, at subsequent times, under a constant- 
magnitude wind blowing from left to right. Since thoroughly burnt areas cannot 
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Figure 2: Basic simulation described by Table O 



be burnt again (on the time scale of the simulation), the area enclosed by the 
front increases with time. The rear, the flanks and the head of the front are 
clearly identifiable. 

The reference simulation is slightly modified to show the ability to deal with 
multiple fronts - Figure [3l It demonstrates the capability to deal with the 
merging of fronts (two main fronts), and to deal with the so-called islands, i.e. 
an unburnt area surrounded by a burnt area. 

In Figures H] and [5l we use the same parameters as in Table [5] but At = 
2.5 • 10~^, and a depends on x, a being equal to 0.5 if a; < 1.7, and a — 0.25 
(Figure 13]) or a = 1.0 (Figure [5]) if x > 1.8, and a being linearly interpolated 
for intermediate values of x. Since a takes into account the available fuel load- 
ing, these two simulations roughly show the influence of the inhomogeneous 
available fuel loading, should it increase (Figured]) or decrease (Figure [5]). The 
inherent decrease of the radius of curvature at the head for a constant-direction 
wind suggests that some vacillation of wind direction contributes when the head 
broadens under otherwise uniform conditions. 

Figure [6] shows the impact of a rotating wind direction. If north is toward 
the top of the figure, then the wind is oriented first west-to-east and tends later 
to south-to-north. 

The next two Figures jT] and [5] show the behavior of two fronts subject to a 
simple-counterflow wind, i.e., a wind defined as: 




(19) 
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0.5 1 1.5 2 2.5 3 

Figure 3: Two main fronts merge, and an island - the unburnt area within the 
biggest front - is burnt. 

2.2 I 1 1 1 1 1 1 1 




0.8 I 1 1 1 1 1 1 1 1 

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

Figure 4: The front slows down at the head for a = 0.25 if a; > 1.8. The final 
time is changed to T/ = 1.5. 
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0.8 I ' ' 1 ' 1 1 

0.5 1 1.5 2 2.5 3 

Figure 5: The front advances faster at the head for a = 1.0 if a; > 1.8. 




0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

Figure 6: Same as the reference simulation, but with a changing wind direction. 
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0.4 

-1.5 -1 -0.5 0.5 1 1.5 

Figure 7: Evolution of the merged front from initially two mirror-image fronts, 
one to each side of the stagnation line for a converging x-component wind, but 
both to one side of the stagnation line for a diverging y-component wind. 

where u is set to 100. A counterflow exemplifies wind conditions well suited for 
setting a backfire, to preburn the vegetation in the path of a wind-aided fire. 

The last Figure [9] shows a front that propagates over an idealized hill. Where 
the slope is positive (between x = 1.6 and x = 1.7), t he firefront typically 
advances faster. Downhill the front typically slows down [Luke and McArthuil . 
19781 pp. 94-97] . The speed function is therefore modified to take into account 
the slope s: 

-^topography = -f^ X 6 , (^0) 

where s is in radians. 

7 Conclusion and Future Prospects 

A semi-empirical, equilibrium-type firespread rate has been used to model a 
wind-aided firefront propagation across wildland surface vegetation. In this for- 
mulation, the rate depends primarily on the wind speed, and the angle between 
the wind direction and the normal to the firefront (idealized as a one-dimensional 
interface). In scenarios arising in practice, the front may consist of several closed 
curves (possibly nested) that can merge as they propagate. 

Level set methods appear capable of treating the model formulated to simu- 
late wildland fire evolution. They treat readily the topological changes that may 
occur to the firefront, and they are known to converge to the physical solution 
of front tracking problem. 
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Figure 8: Evolution of the merged front from initially two mirror-image fronts, 
here symmetrically sited relative to a simple counterflow wind. 
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0.8 I 1 1 1 1 1 1 1 1 

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

Figure 9: Taking into account topography: the front propagates over an ideal- 
ized hill. 

They were applied via the Multivac package. This open-source library is 
designed to handle a wide range of applications without loss of computing per- 
formance. It includes several algorithms and numerical schemes, primarily for 
the narrow-band level set method, which is more computationally efficient than 
the full level set method. 

A possible direction for future work is to focus on parameter estimation 
within the context of the simple model illustrated herein. A cost is introduced 
to measure the distance between the simulated front and ground, aerial, and/or 
satellite observations. The discrepancy between the simulated and observed po- 
sitions of the front may be based either on the front arrival times (at monitored 
locations), or on distances between the simulated front and the monitored lo- 
cations (at arrival times). For gradient-based optimization methods, the main 
challenge is to compute the derivative of the cost function with respect to the 
parameters. An adjoint code being difficult to construct, alternative methods 
should be sought. 

This work could help guide fire-control tactics. The objective function would 
then penalize front advance into societal assets, and penalize the cost of the 
firefighting activity. The parameters would be the model variables modifiable 
by firefighting countermeasures. The links between this optimization problem 
and shape optimization should be investigated. 
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